function [err, FPresult]=FP_TC(DFILE, sf, plot_flag)
% Analyze Field Potential time course

global VOLTAGE CURRENT
global CONTROL

%
err = 0;

[DFILE] = synch_file(DFILE, sf);
if(isempty(DFILE) | isempty(VOLTAGE))
   fprintf('? could not get dfile');
   return;
end;

QueMessage('FP_TC - Starting');
dat = [];
time = [];
ivresult = []; % initialize it
do_spike = 0;	% 0 turns off detailed spike analysis; 1 turns it on 

update = 0; % flag
if(~check_field(CONTROL(sf),'stim_time')) % doesn't exist - then create it.
   CONTROL(sf).stim_time=[5 55];
   update = 1;
end;
if(~check_field(CONTROL(sf),'psp_time'))
   CONTROL(sf).psp_time= [9 11 31, 9 51 71];
   update=1;
end
if(~check_field(CONTROL(sf),'psp_reclist'))
   CONTROL(sf).psp_reclist= '';
   update=1;
end
if(update==1)
   QueMessage('Updated Control with psp params');
end
% abstract general information
QueMessage('Preparing for analysis');
protocol=deblank(lower(CONTROL(sf).protocol));
rate = DFILE.rate*DFILE.nr_channel/1000;

[records,pts]=size(VOLTAGE);
% always compute the time of each stimulus, in seconds from the start of the data
a = DFILE.ztime;
TM = (a - a(1))/1000; % time at which data was collected... gathered from ztime...
if(max(TM) == 0)
   TM = 15*[0:records-1];
   rate = rate*2;
end;

% Get analysis windows 
if(ischar(CONTROL(sf).stim_time))
   pspt = CONTROL(sf).stim_time;
   p = find( pspt~= ' ');
   pspt = pspt(p);
   stim_time=seq_parse(pspt);; % get the stimulus array - usually 2 times....
   stim_time = stim_time{:};
else
   stim_time = CONTROL(sf).stim_time;
end;

% get the psp definition array. For field potentials, we have the following
% sets of times: (layer 2 DCN pf FP)
% It is assumed that the same response is repeated with similar latencies for each stimulus...

if(ischar(CONTROL(sf).psp_time))
   pspt = CONTROL(sf).psp_time;
   p = find(pspt~= ' ');
   pspt = pspt(p);
   psp_time=seq_parse(pspt); % get psp definition array : for 
   psp_time = psp_time{:};
else
   psp_time = CONTROL(sf).psp_time;
end;


ts1=number_arg(stim_time(1)); % get the baseline duration (most of time before stimulus)
if(length(number_arg(stim_time) > 1))
   ts2=number_arg(stim_time(2)); % get the baseline duration (most of time before stimulus)
   stim_intvl = ts2-ts1;
else
   stim_intvl = 0;
end;

tn1b=number_arg(psp_time(1)); % get psp list - n1 window (measure minimum)
tn1e=number_arg(psp_time(2));

tn2b=number_arg(psp_time(3)); % get psp list - n2 window (measure minimum)
tn2e=number_arg(psp_time(4));

tp3b=number_arg(psp_time(5)); % get psp list - p3 window (measure maximum)
tp3e=number_arg(psp_time(6));

if(length(psp_time) >= 8)
   tp4b=number_arg(psp_time(7)); % get psp list - p4 window (measure minimum)
   tp4e=number_arg(psp_time(8));
else
   tp4b = 0;
   tp4e = 0;
end;
if(stim_intvl > 0)
   tn1b2=tn1b+stim_intvl;
   tn1e2=tn1e+stim_intvl;
   tn2b2=tn2b+stim_intvl;
   tn2e2=tn2e+stim_intvl;
   tp3b2=tp3b+stim_intvl;
   tp3e2=tp3e+stim_intvl;
   if(length(psp_time) >= 8)
      tp4b2=tp4b+stim_intvl;
      tp4e2=tp4e+stim_intvl;
   end;
end;
% compute the time base for plotting (time is [rec, npts] array for EACH record)
time=make_time(DFILE);
tmax=max(max(time));
k = find(DFILE.rate < 1);
if(DFILE.mode <= 2)
   DFILE.rate(k) = 100;
end;
RATES = (DFILE.rate * DFILE.nr_channel)/ 1000; % array of sampling rates, convert to msec
if(DFILE.mode <= 2)
   RATES=ones(length(RATES), 1)*0.2;
end;

% calculate measurement times
for i=1:records
   fp_base(i) = floor((ts1/RATES(i) - 0.5)*0.75); % for RMP/ihold determination
   
   fp_n1b(i)=floor(tn1b/RATES(i)); % for peak n1
   fp_n1e(i)=floor(tn1e/RATES(i)); % for peak n1
   
   fp_n2b(i)=floor(tn2b/RATES(i)); % for peak n2
   fp_n2e(i)=floor(tn2e/RATES(i)); % for peak n2
   
   fp_p3b(i)=floor(tp3b/RATES(i)); % for peak p3
   fp_p3e(i)=floor(tp3e/RATES(i)); % for peak p3
   
   if(length(psp_time) >= 8)
      fp_p4b(i)=floor(tp4b/RATES(i)); % for peak p4
      fp_p4e(i)=floor(tp4e/RATES(i)); % for peak p4	
   end;
   
   if(stim_intvl > 0)
      fp_n1b2(i)=floor(tn1b2/RATES(i)); % for peak n1
      fp_n1e2(i)=floor(tn1e2/RATES(i)); % for peak n1
      
      fp_n2b2(i)=floor(tn2b2/RATES(i)); % for peak n2
      fp_n2e2(i)=floor(tn2e2/RATES(i)); % for peak n2
      
      fp_p3b2(i)=floor(tp3b2/RATES(i)); % for peak p3
      fp_p3e2(i)=floor(tp3e2/RATES(i)); % for peak p3
      
      if(length(psp_time) >= 8)
         fp_p4b2(i)=floor(tp4b2/RATES(i)); % for peak p4
         fp_p4e2(i)=floor(tp4e2/RATES(i)); % for peak p4
      end;
   end;
end;

% measure holding current and "rmp"
for i = 1:records
   hold_cur(i) = mean(CURRENT(i,1:fp_base(i))); % save in array for later usage
   rmp_volt(i) = mean(VOLTAGE(i,1:fp_base(i)));
end;
CONTROL(sf).iHold = mean(hold_cur);
CONTROL(sf).Rmp = mean(rmp_volt);

% but first smooth the voltage out a bit
for i = 1:records
   fsamp = 1000/RATES(i); % get sampling frequency
   fco = 500;		% cutoff frequency in Hz
   wco = fco/(fsamp/2); % wco of 1 is for half of the sample rate, so set it like this...
   if(wco < 1) % if wco is > 1 then this is not a filter!
      [b, a] = fir_win(8, wco); % fir type filter... seems to work best, with highest order min distortion of dv/dt...
      vsmo(i,:) = DigitalFilt(b, a, VOLTAGE(i,:)')'; % filter all the traces...
   else
      vsmo(i,:) = VOLTAGE(i,:);
   end
end

for i = 1:records
   base = mean(vsmo(i,1:fp_base(i)));
   vsmo(i,:) = vsmo(i,:) - base;
end;
QueMessage('FP_TC - Measureing FP peak voltages and P3 slope for stim 1');
% find min voltage in window 1 (N2, pop spike)
MAX_N1V = zeros(records, 1);
MAX_N2V = zeros(records, 1);
MAX_P3V = zeros(records, 1);
MAX_P3S = zeros(records, 1);

if(stim_intvl > 0)
   MAX_N1V2 = zeros(records, 1);
   MAX_N2V2 = zeros(records, 1);
   MAX_P3V2 = zeros(records, 1);
   MAX_P3S2 = zeros(records, 1);
end;

for i = 1:records
   bdpts = any(abs(rmp_volt(i)) > 15); % 1 mv baseline? bad points!
   MAX_N1V(i) = -min(vsmo(i,fp_n1b(i):fp_n1e(i))'); % voltage peak N1, 	flipped for direct comparison...
   MAX_N1V(i,bdpts) = NaN;
   
   MAX_N2V(i) = -min(vsmo(i,fp_n2b(i):fp_n2e(i))'); % voltage peak N2
   MAX_N2V(i,bdpts) = NaN;
   
   [MAX_P3V(i), pk] = max(vsmo(i, fp_p3b(i):fp_p3e(i))'); % voltage peak P3
   MAX_P3V(i,bdpts) = NaN;
   
   MAX_P3S(i) = max(diff(vsmo(i,fp_p3b(i):(fp_p3b(i)+pk))')) / RATES(i); % max slope on rising phase
   MAX_P3S(i,bdpts) = NaN;
   if(stim_intvl > 0)
      MAX_N1V2(i) = -min(vsmo(i,fp_n1b2(i):fp_n1e2(i))'); % voltage peak N1, 	flipped for direct comparison...
      MAX_N1V2(i,bdpts) = NaN;
      
      MAX_N2V2(i) = -min(vsmo(i,fp_n2b2(i):fp_n2e2(i))'); % voltage peak N2
      MAX_N2V2(i,bdpts) = NaN;
      
      [MAX_P3V2(i), pk] = max(vsmo(i, fp_p3b2(i):fp_p3e2(i))'); % voltage peak P3
      MAX_P3V2(i,bdpts) = NaN;
      
      MAX_P3S2(i) = max(diff(vsmo(i,fp_p3b2(i):(fp_p3b2(i)+pk))')) / RATES(i); % max slope on rising phase
      MAX_P3S2(i,bdpts) = NaN;
   end;
end;

FP.maxN1=MAX_N1V;
FP.maxN2=MAX_N2V;
FP.maxP3=MAX_P3V;
FP.maxP3s=MAX_P3S;
FP.maxN12=MAX_N1V2;
FP.maxN22=MAX_N2V2;
FP.maxP32=MAX_P3V2;
FP.maxP3s2=MAX_P3S2;
FP.time = TM;
% do F test on the windows...
%tbs=1;
%a=find(TM<600);
%tbe=a(length(a)); % use the last point accessible
%a=find(TM>900);
%tps=a(1);
%a=find(TM>1500);
%tpe=a(1);
%tbs
%tbe
%tps
%tpe
%size(TM)

%[P, F, dfR, dfF, yp1, yp2, xf, yf] = bears_f_test(TM(tbs:tbe)', MAX_FP2V(tbs:tbe), ...
%   TM(tps:tpe)', MAX_FP2V(tps:tpe));
%disp(sprintf('base: %7.2f mV  Sig: %7.2f mV, P=%6.3f', mean(MAX_FP2V(tbs:tbe)), ...
%   mean(MAX_FP2V(tps:tpe)), P))

CONTROL(sf).FP = FP;

QueMessage('EPSP_TC - analysis complete');

% plot if figure is set

h = findobj('Tag', 'FP_TC'); % check for pre-existing window
if(isempty(h)) % if none, make one
   h = figure('Tag', 'FP_TC', 'Name', 'FP Time Course Analysis', 'NumberTitle', 'off');
end
figure(h); % otherwise, select it
clf; % always clear the window...
fsize = 7;
msize = 3;
tmax = max(time);    

%plot Min voltage for N2
subplot('Position',[0.07,0.70,0.7,0.250]);
plot(TM(1:records),MAX_N1V(1:records), '-ks', 'Markersize', msize, 'Markerfacecolor', 'k') % data in black
hold on;
plot(TM(1:records), MAX_N2V(1:records), '-b^', 'Markersize', msize, 'Markerfacecolor', 'b'); % data in blue
plot(TM(1:records), MAX_P3V(1:records), '-ro', 'Markersize', msize, 'Markerfacecolor', 'r'); % data in red

% open symbols for second psp.
if(stim_intvl > 0)
   plot(TM(1:records),MAX_N1V2(1:records), '-ks', 'Markersize', msize, 'Markerfacecolor', 'w') % data in black
   plot(TM(1:records), MAX_N2V2(1:records), '-b^', 'Markersize', msize, 'Markerfacecolor', 'w'); % data in blue
   plot(TM(1:records), MAX_P3V2(1:records), '-ro', 'Markersize', msize, 'Markerfacecolor', 'w'); % data in red
end;

set(gca, 'FontSize', fsize);
ylabel('FP Amplitude (mV)');
grid;
u=get(gca, 'YLim');
hold off;
%plot Max Slope voltage 
subplot('Position',[0.07,0.40,0.7,0.250]);
plot(TM(1:records),MAX_P3S(1:records), '-r^', 'Markersize', msize, 'Markerfacecolor', 'r') % data in red
set(gca, 'FontSize', fsize);
v=get(gca, 'YLim');
set(gca, 'Ylim', [0 v(2)]);
grid;
ylabel('FP P3 rising Slope (mV/ms)');

subplot('Position', [0.07, 0.1, 0.7, 0.25]);
plot(TM(1:records), MAX_N1V2./MAX_N1V, '-ks', 'Markersize', msize, 'Markerfacecolor', 'k');
plot(TM(1:records), MAX_N2V2./MAX_N2V, '-bx', 'Markersize', msize, 'Markerfacecolor', 'b');
plot(TM(1:records), MAX_P3V2./MAX_P3V, '-ro', 'Markersize', msize, 'Markerfacecolor', 'r');

% control printing and closing of window for automatic runs
% f = 1 creates plot and leaves it up
% f = 2 creates plot but closes it when done
% f = 3 creates plot and prints it and then closes it
% control printing and closing of window for automatic runs
% plot_flag = 0 creates plot and leaves it up
% plot_flag = 1 creates plot and prints to the printer
% plot_flag = 2 creates plot and prints it and then closes it
% checking the pflag box in datac main window is plot_flag = 1;
if (plot_flag > 0)
   h = findobj('Tag', 'FP_IO');
   hg = findobj('Tag', 'FP_info');
   set(hg, 'Visible', 'off');
   print(h);
   set(hg, 'Visible', 'on');
end
if plot_flag == 2  
   close
end


